To-do

Rather than focus on p-values, do something like a confusion matrix. 1) set up a cutoff of loadings and/or VIP scores that make a variable a significant contributor. 2) Quantify how often each method gets this right (was it a discriminating variable vs. was it detected as a discriminating variable) 3) summarize this as proportion of false positives or some other index (Ask Matt, or look at notes) # Introduction The purpose of this notebook is to stress-test PCA, RDA, and PLS-DA to better understand their properties. All the datasets in this notebook have 35 variables and 20 observations split into two groups (“a” and “b”). I randomly generate nperm number of data frames using the same parameters, then obtain p-values from each statistical method on all data sets to test if the groups differ. I then visualize the distribution of p-values.

Packages needed

library(chemhelper) # for sim_*() functions
library(tidyverse)
library(ropls) # I'm open to using pls or another package to fit PLS-DA models.
library(iheatmapr) # for interactive correlation heatmaps
# library(broom)
library(vegan) # for rda()

Define Functions

PCA t-test

Does PCA, then does t-test on PC1 scores, returns p-value.

PCA.t.test <- function(df){
  pca.m <- opls(select(df, -group), plotL = FALSE, printL = FALSE)
  pca.d <- get_plotdata(pca.m)$plot_data
  pca.t <- t.test(pca.d$p1 ~ df$group)
  return(pca.t$p.value)
}

RDA ANOVA

Fits RDA model, runs ANOVA for global significance of model, returns p-value.

RDA.anova <- function(df){
  rda.m <- rda(select(df, -group) ~ df$group)
  return(anova(rda.m)$`Pr(>F)`[[1]])
}

PLS-DA

Fits PLS-DA model, calculates p-values with 100 permutations, returns \(Q^2\), \(p_{(R^2_Y)}\), and \(p_{(Q^2)}\)

PLS.test <- function(df){
  pls.m <- opls(select(df, -group), df$group, permI = 100, predI = 2, printL = FALSE, plotL = FALSE)
  return(select(pls.m@summaryDF, `Q2(cum)`, pR2Y, pQ2))
}

1. Many, highly multicollinear distinguishing variables contributing to large effect size.

This code creates a dataset with 30 discriminating variables that co-vary and 5 uncorrelated variables.

base_df <- sim_df(20, 2)
nperm = 50
df1 <- map(1:nperm,
    ~base_df %>% 
      sim_covar(5, var = 1, cov = 0, name = "uncorr") %>% 
      sim_discr(30, var = 1, cov = 0.5, group_means = c(-0.5, 0.5), name = "discr")
) %>% set_names(paste("dataset", 1:nperm))
df1[[1]] %>% 
  select(-group) %>% 
  cor() %>% 
  iheatmap(row_labels = TRUE, col_labels = TRUE, name = "cor")

PCA 1

pvals.pca1 <- map_dbl(df1, PCA.t.test)

RDA 1

pvals.rda1 <- map_dbl(df1, RDA.anova)

PLS 1

output.pls1 <- map(df1, PLS.test) %>% bind_rows()

Results 1

results1 <- output.pls1 %>% add_column(pvals.rda1, pvals.pca1)
results1 %>% summarize_all(funs(mean = mean))
pvals1 <- results1 %>% 
  select(`PLS pQ2` = pQ2, `RDA ANOVA` = pvals.rda1, `PCA t-test` = pvals.pca1) %>% 
  gather(key = "Method", value = "p value")
ggplot(pvals1, aes(x = `p value`, fill = Method)) +
  # geom_density(alpha = 0.3) +
  geom_histogram(bins = 20, position = "dodge") +
  geom_vline(aes(xintercept = 0.05), color = "red")

In this scenario, PLS gives a low p-value. Since all 30 of the discriminating variables covary strongly with eachother, I interpret this as PLS being highly robust to multicollinearity. That is, there aren’t really 30 discriminating variables, but one latent variable that affects all 30 variables. I imagine that if I were to create 30 discriminating variables in 6 groups of 5 variables that the p-values of all methods would be similar.

2. Needle in a haystack—many, highly collinear variables, few discriminating variables

df2 <- map(1:nperm, ~base_df %>% 
  sim_covar(10, 1, 0.5, "corrA") %>% 
  sim_covar(10, 1, 0.5, "corrB") %>% 
  sim_covar(10, 1, 0.5, "corrC") %>% 
  sim_discr(5, 1, 0, group_means = c(-1, 1), name = "discr")
  )
df2[[2]] %>% 
  select(-group) %>% 
  cor() %>% 
  iheatmap(row_labels = TRUE, col_labels = TRUE, name = "cor")

PCA 2

pvals.pca2 <- map_dbl(df2, PCA.t.test)
mean(pvals.pca2)
[1] 0.3589526

RDA 2

pvals.rda2 <- map_dbl(df2, RDA.anova)
mean(pvals.rda2)
[1] 0.0019

PLS 2

output.pls2 <- map(df2, PLS.test) %>% bind_rows()

Results 2

results2 <- output.pls2 %>% add_column(pvals.rda2, pvals.pca2)
results2 %>% summarize_all(funs(mean = mean))
pvals2 <- results2 %>% 
  select(`PLS pQ2` = pQ2, `RDA ANOVA` = pvals.rda2, `PCA t-test` = pvals.pca2) %>% 
  gather(key = "Method", value = "p value")
ggplot(pvals2, aes(x = `p value`, fill = Method)) +
  # geom_density(alpha = 0.3) + 
  geom_histogram(position = "dodge", bins = 20) +
  # ylim(0, 20) +
  # xlim(0, 0.5) +
  geom_vline(aes(xintercept = 0.05), color = "red")

Here PCA does less well because the variables that contribute most to variation in the dataset are not those variables that contribute to differences between groups. PLS-DA and RDA both consistently find differences between the groups. However, these significant differences are only due to 5 variables. It’s important to note how much variation is explained by the predictive axes of PLS or the constrained portion of RDA to interpret these results.

3. Red Herring? (AKA Needle in a highly multicollinear haystack.)

Ok, this is ridiculous, but here is a dataset with 30, strongly multicollinear variables that don’t discriminate between groups and only 5 variables that do discriminate between groups.

df3 <- map(1:nperm,
             ~base_df %>% 
               sim_covar(30, 1, 0.7, name = "corr") %>% 
               sim_discr(5, 1, 0, name = "discr", group_means = c(-0.5, 0.5))
)
df3[[1]] %>% 
  select(-group) %>% 
  cor() %>% 
  iheatmap(row_labels = TRUE, col_labels = TRUE, name = "cor")

PCA 3

pvals.pca3 <- map_dbl(df3, PCA.t.test)
mean(pvals.pca3)
[1] 0.4606713

RDA 2

pvals.rda3 <- map_dbl(df3, RDA.anova)
mean(pvals.rda3)
[1] 0.19278

PLS 2

output.pls3 <- map(df3, PLS.test) %>% bind_rows()

Results 3

results3 <- output.pls3 %>% add_column(pvals.rda3, pvals.pca3)
results3 %>% summarize_all(funs(mean = mean))
meanQ <- mean(results3$`Q2(cum)`)
pvals3 <- results3 %>% 
  select(`PLS pQ2` = pQ2, `RDA ANOVA` = pvals.rda3, `PCA t-test` = pvals.pca3) %>% 
  gather(key = "Method", value = "p value")
ggplot(pvals3, aes(x = `p value`, fill = Method)) +
  # geom_density(alpha = 0.3) + 
  geom_histogram(position = "dodge", bins = 20) +
  geom_vline(aes(xintercept = 0.05), color = "red")

Even in this ridiculous scenario, PLS is able to “find” the 5 discriminating variables most of the time and give a significant p-value. The mean \(Q^2\) value is 0.39, so it’s possible that some of those “significant” p-values were obtained from models with poor predictive power. Because within group co-variation is high and between-group co-variation is low, the F-test used to analyze RDA results gives a high p-value. I think.

LS0tCnRpdGxlOiAiU3RyZXNzIHRlc3RpbmcgUENBLCBSREEsIGFuZCBQTFMiCmF1dGhvcjogIkVyaWMgUi4gU2NvdHQiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazogCiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKICAgIHRvYzogeWVzCi0tLQoKIyBUby1kbwpSYXRoZXIgdGhhbiBmb2N1cyBvbiBwLXZhbHVlcywgZG8gc29tZXRoaW5nIGxpa2UgYSBjb25mdXNpb24gbWF0cml4LiAgMSkgc2V0IHVwIGEgY3V0b2ZmIG9mIGxvYWRpbmdzIGFuZC9vciBWSVAgc2NvcmVzIHRoYXQgbWFrZSBhIHZhcmlhYmxlIGEgc2lnbmlmaWNhbnQgY29udHJpYnV0b3IuICAyKSBRdWFudGlmeSBob3cgb2Z0ZW4gZWFjaCBtZXRob2QgZ2V0cyB0aGlzIHJpZ2h0ICh3YXMgaXQgYSBkaXNjcmltaW5hdGluZyB2YXJpYWJsZSB2cy4gd2FzIGl0IGRldGVjdGVkIGFzIGEgZGlzY3JpbWluYXRpbmcgdmFyaWFibGUpIDMpIHN1bW1hcml6ZSB0aGlzIGFzIHByb3BvcnRpb24gb2YgZmFsc2UgcG9zaXRpdmVzIG9yIHNvbWUgb3RoZXIgaW5kZXggKEFzayBNYXR0LCBvciBsb29rIGF0IG5vdGVzKQojIEludHJvZHVjdGlvbgpUaGUgcHVycG9zZSBvZiB0aGlzIG5vdGVib29rIGlzIHRvIHN0cmVzcy10ZXN0IFBDQSwgUkRBLCBhbmQgUExTLURBIHRvIGJldHRlciB1bmRlcnN0YW5kIHRoZWlyIHByb3BlcnRpZXMuICBBbGwgdGhlIGRhdGFzZXRzIGluIHRoaXMgbm90ZWJvb2sgaGF2ZSAzNSB2YXJpYWJsZXMgYW5kIDIwIG9ic2VydmF0aW9ucyBzcGxpdCBpbnRvIHR3byBncm91cHMgKCJhIiBhbmQgImIiKS4gSSByYW5kb21seSBnZW5lcmF0ZSBgbnBlcm1gIG51bWJlciBvZiBkYXRhIGZyYW1lcyB1c2luZyB0aGUgc2FtZSBwYXJhbWV0ZXJzLCB0aGVuIG9idGFpbiBwLXZhbHVlcyBmcm9tIGVhY2ggc3RhdGlzdGljYWwgbWV0aG9kIG9uIGFsbCBkYXRhIHNldHMgdG8gdGVzdCBpZiB0aGUgZ3JvdXBzIGRpZmZlci4gIEkgdGhlbiB2aXN1YWxpemUgdGhlIGRpc3RyaWJ1dGlvbiBvZiBwLXZhbHVlcy4KCiMgUGFja2FnZXMgbmVlZGVkCgpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGNoZW1oZWxwZXIpICMgZm9yIHNpbV8qKCkgZnVuY3Rpb25zCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHJvcGxzKSAjIEknbSBvcGVuIHRvIHVzaW5nIHBscyBvciBhbm90aGVyIHBhY2thZ2UgdG8gZml0IFBMUy1EQSBtb2RlbHMuCmxpYnJhcnkoaWhlYXRtYXByKSAjIGZvciBpbnRlcmFjdGl2ZSBjb3JyZWxhdGlvbiBoZWF0bWFwcwojIGxpYnJhcnkoYnJvb20pCmxpYnJhcnkodmVnYW4pICMgZm9yIHJkYSgpCmBgYAoKIyBEZWZpbmUgRnVuY3Rpb25zCgojIyBQQ0EgdC10ZXN0CgpEb2VzIFBDQSwgdGhlbiBkb2VzIHQtdGVzdCBvbiBQQzEgc2NvcmVzLCByZXR1cm5zIHAtdmFsdWUuCgpgYGB7cn0KUENBLnQudGVzdCA8LSBmdW5jdGlvbihkZil7CiAgcGNhLm0gPC0gb3BscyhzZWxlY3QoZGYsIC1ncm91cCksIHBsb3RMID0gRkFMU0UsIHByaW50TCA9IEZBTFNFKQogIHBjYS5kIDwtIGdldF9wbG90ZGF0YShwY2EubSkkcGxvdF9kYXRhCiAgcGNhLnQgPC0gdC50ZXN0KHBjYS5kJHAxIH4gZGYkZ3JvdXApCiAgcmV0dXJuKHBjYS50JHAudmFsdWUpCn0KYGBgCgojIyBSREEgQU5PVkEKCkZpdHMgUkRBIG1vZGVsLCBydW5zIEFOT1ZBIGZvciBnbG9iYWwgc2lnbmlmaWNhbmNlIG9mIG1vZGVsLCByZXR1cm5zIHAtdmFsdWUuCgpgYGB7cn0KUkRBLmFub3ZhIDwtIGZ1bmN0aW9uKGRmKXsKICByZGEubSA8LSByZGEoc2VsZWN0KGRmLCAtZ3JvdXApIH4gZGYkZ3JvdXApCiAgcmV0dXJuKGFub3ZhKHJkYS5tKSRgUHIoPkYpYFtbMV1dKQp9CmBgYAoKIyMgUExTLURBCgpGaXRzIFBMUy1EQSBtb2RlbCwgY2FsY3VsYXRlcyBwLXZhbHVlcyB3aXRoIDEwMCBwZXJtdXRhdGlvbnMsIHJldHVybnMgJFFeMiQsICRwX3soUl4yX1kpfSQsIGFuZCAkcF97KFFeMil9JAoKYGBge3J9ClBMUy50ZXN0IDwtIGZ1bmN0aW9uKGRmKXsKICBwbHMubSA8LSBvcGxzKHNlbGVjdChkZiwgLWdyb3VwKSwgZGYkZ3JvdXAsIHBlcm1JID0gMTAwLCBwcmVkSSA9IDIsIHByaW50TCA9IEZBTFNFLCBwbG90TCA9IEZBTFNFKQogIHJldHVybihzZWxlY3QocGxzLm1Ac3VtbWFyeURGLCBgUTIoY3VtKWAsIHBSMlksIHBRMikpCn0KYGBgCgoKCiMgMS4gTWFueSwgaGlnaGx5IG11bHRpY29sbGluZWFyIGRpc3Rpbmd1aXNoaW5nIHZhcmlhYmxlcyBjb250cmlidXRpbmcgdG8gbGFyZ2UgZWZmZWN0IHNpemUuCgpUaGlzIGNvZGUgY3JlYXRlcyBhIGRhdGFzZXQgd2l0aCAzMCBkaXNjcmltaW5hdGluZyB2YXJpYWJsZXMgdGhhdCBjby12YXJ5IGFuZCA1IHVuY29ycmVsYXRlZCB2YXJpYWJsZXMuCgpgYGB7cn0KYmFzZV9kZiA8LSBzaW1fZGYoMjAsIDIpCm5wZXJtID0gNTAKYGBgCgpgYGB7cn0KZGYxIDwtIG1hcCgxOm5wZXJtLAogICAgfmJhc2VfZGYgJT4lIAogICAgICBzaW1fY292YXIoNSwgdmFyID0gMSwgY292ID0gMCwgbmFtZSA9ICJ1bmNvcnIiKSAlPiUgCiAgICAgIHNpbV9kaXNjcigzMCwgdmFyID0gMSwgY292ID0gMC41LCBncm91cF9tZWFucyA9IGMoLTAuNSwgMC41KSwgbmFtZSA9ICJkaXNjciIpCikgJT4lIHNldF9uYW1lcyhwYXN0ZSgiZGF0YXNldCIsIDE6bnBlcm0pKQoKZGYxW1sxXV0gJT4lIAogIHNlbGVjdCgtZ3JvdXApICU+JSAKICBjb3IoKSAlPiUgCiAgaWhlYXRtYXAocm93X2xhYmVscyA9IFRSVUUsIGNvbF9sYWJlbHMgPSBUUlVFLCBuYW1lID0gImNvciIpCmBgYAoKKipQQ0EgMSoqCgpgYGB7cn0KcHZhbHMucGNhMSA8LSBtYXBfZGJsKGRmMSwgUENBLnQudGVzdCkKYGBgCgoqKlJEQSAxKioKCmBgYHtyfQpwdmFscy5yZGExIDwtIG1hcF9kYmwoZGYxLCBSREEuYW5vdmEpCmBgYAoKKipQTFMgMSoqCgpgYGB7cn0Kb3V0cHV0LnBsczEgPC0gbWFwKGRmMSwgUExTLnRlc3QpICU+JSBiaW5kX3Jvd3MoKQpgYGAKCiMjIFJlc3VsdHMgMQoKYGBge3J9CnJlc3VsdHMxIDwtIG91dHB1dC5wbHMxICU+JSBhZGRfY29sdW1uKHB2YWxzLnJkYTEsIHB2YWxzLnBjYTEpCnJlc3VsdHMxICU+JSBzdW1tYXJpemVfYWxsKGZ1bnMobWVhbiA9IG1lYW4pKQpgYGAKCmBgYHtyfQpwdmFsczEgPC0gcmVzdWx0czEgJT4lIAogIHNlbGVjdChgUExTIHBRMmAgPSBwUTIsIGBSREEgQU5PVkFgID0gcHZhbHMucmRhMSwgYFBDQSB0LXRlc3RgID0gcHZhbHMucGNhMSkgJT4lIAogIGdhdGhlcihrZXkgPSAiTWV0aG9kIiwgdmFsdWUgPSAicCB2YWx1ZSIpCgpnZ3Bsb3QocHZhbHMxLCBhZXMoeCA9IGBwIHZhbHVlYCwgZmlsbCA9IE1ldGhvZCkpICsKICAjIGdlb21fZGVuc2l0eShhbHBoYSA9IDAuMykgKwogIGdlb21faGlzdG9ncmFtKGJpbnMgPSAyMCwgcG9zaXRpb24gPSAiZG9kZ2UiKSArCiAgZ2VvbV92bGluZShhZXMoeGludGVyY2VwdCA9IDAuMDUpLCBjb2xvciA9ICJyZWQiKQoKYGBgCgpJbiB0aGlzIHNjZW5hcmlvLCBQTFMgZ2l2ZXMgYSBsb3cgcC12YWx1ZS4gIFNpbmNlIGFsbCAzMCBvZiB0aGUgZGlzY3JpbWluYXRpbmcgdmFyaWFibGVzIGNvdmFyeSBzdHJvbmdseSB3aXRoIGVhY2hvdGhlciwgSSBpbnRlcnByZXQgdGhpcyBhcyBQTFMgYmVpbmcgaGlnaGx5IHJvYnVzdCB0byBtdWx0aWNvbGxpbmVhcml0eS4gIFRoYXQgaXMsIHRoZXJlIGFyZW4ndCAqcmVhbGx5KiAzMCBkaXNjcmltaW5hdGluZyB2YXJpYWJsZXMsIGJ1dCBvbmUgbGF0ZW50IHZhcmlhYmxlIHRoYXQgYWZmZWN0cyBhbGwgMzAgdmFyaWFibGVzLiAgSSBpbWFnaW5lIHRoYXQgaWYgSSB3ZXJlIHRvIGNyZWF0ZSAzMCBkaXNjcmltaW5hdGluZyB2YXJpYWJsZXMgaW4gNiBncm91cHMgb2YgNSB2YXJpYWJsZXMgdGhhdCB0aGUgcC12YWx1ZXMgb2YgYWxsIG1ldGhvZHMgd291bGQgYmUgc2ltaWxhci4KCiMgMi4gTmVlZGxlIGluIGEgaGF5c3RhY2stLS1tYW55LCBoaWdobHkgY29sbGluZWFyIHZhcmlhYmxlcywgZmV3IGRpc2NyaW1pbmF0aW5nIHZhcmlhYmxlcwoKYGBge3J9CmRmMiA8LSBtYXAoMTpucGVybSwgfmJhc2VfZGYgJT4lIAogIHNpbV9jb3ZhcigxMCwgMSwgMC41LCAiY29yckEiKSAlPiUgCiAgc2ltX2NvdmFyKDEwLCAxLCAwLjUsICJjb3JyQiIpICU+JSAKICBzaW1fY292YXIoMTAsIDEsIDAuNSwgImNvcnJDIikgJT4lIAogIHNpbV9kaXNjcig1LCAxLCAwLCBncm91cF9tZWFucyA9IGMoLTEsIDEpLCBuYW1lID0gImRpc2NyIikKICApCgpkZjJbWzJdXSAlPiUgCiAgc2VsZWN0KC1ncm91cCkgJT4lIAogIGNvcigpICU+JSAKICBpaGVhdG1hcChyb3dfbGFiZWxzID0gVFJVRSwgY29sX2xhYmVscyA9IFRSVUUsIG5hbWUgPSAiY29yIikKYGBgCgoqKlBDQSAyKioKCmBgYHtyfQpwdmFscy5wY2EyIDwtIG1hcF9kYmwoZGYyLCBQQ0EudC50ZXN0KQptZWFuKHB2YWxzLnBjYTIpCmBgYAoKKipSREEgMioqCgpgYGB7cn0KcHZhbHMucmRhMiA8LSBtYXBfZGJsKGRmMiwgUkRBLmFub3ZhKQptZWFuKHB2YWxzLnJkYTIpCmBgYAoKKipQTFMgMioqCgpgYGB7cn0Kb3V0cHV0LnBsczIgPC0gbWFwKGRmMiwgUExTLnRlc3QpICU+JSBiaW5kX3Jvd3MoKQpgYGAKCiMjIFJlc3VsdHMgMgoKYGBge3J9CnJlc3VsdHMyIDwtIG91dHB1dC5wbHMyICU+JSBhZGRfY29sdW1uKHB2YWxzLnJkYTIsIHB2YWxzLnBjYTIpCnJlc3VsdHMyICU+JSBzdW1tYXJpemVfYWxsKGZ1bnMobWVhbiA9IG1lYW4pKQpgYGAKCmBgYHtyfQpwdmFsczIgPC0gcmVzdWx0czIgJT4lIAogIHNlbGVjdChgUExTIHBRMmAgPSBwUTIsIGBSREEgQU5PVkFgID0gcHZhbHMucmRhMiwgYFBDQSB0LXRlc3RgID0gcHZhbHMucGNhMikgJT4lIAogIGdhdGhlcihrZXkgPSAiTWV0aG9kIiwgdmFsdWUgPSAicCB2YWx1ZSIpCgpnZ3Bsb3QocHZhbHMyLCBhZXMoeCA9IGBwIHZhbHVlYCwgZmlsbCA9IE1ldGhvZCkpICsKICAjIGdlb21fZGVuc2l0eShhbHBoYSA9IDAuMykgKyAKICBnZW9tX2hpc3RvZ3JhbShwb3NpdGlvbiA9ICJkb2RnZSIsIGJpbnMgPSAyMCkgKwogICMgeWxpbSgwLCAyMCkgKwogICMgeGxpbSgwLCAwLjUpICsKICBnZW9tX3ZsaW5lKGFlcyh4aW50ZXJjZXB0ID0gMC4wNSksIGNvbG9yID0gInJlZCIpCmBgYApIZXJlIFBDQSBkb2VzIGxlc3Mgd2VsbCBiZWNhdXNlIHRoZSB2YXJpYWJsZXMgdGhhdCBjb250cmlidXRlIG1vc3QgdG8gdmFyaWF0aW9uIGluIHRoZSBkYXRhc2V0IGFyZSBub3QgdGhvc2UgdmFyaWFibGVzIHRoYXQgY29udHJpYnV0ZSB0byBkaWZmZXJlbmNlcyBiZXR3ZWVuIGdyb3Vwcy4gUExTLURBIGFuZCBSREEgYm90aCBjb25zaXN0ZW50bHkgZmluZCBkaWZmZXJlbmNlcyBiZXR3ZWVuIHRoZSBncm91cHMuICBIb3dldmVyLCB0aGVzZSBzaWduaWZpY2FudCBkaWZmZXJlbmNlcyBhcmUgb25seSBkdWUgdG8gKio1KiogdmFyaWFibGVzLiAgSXQncyBpbXBvcnRhbnQgdG8gbm90ZSBob3cgbXVjaCB2YXJpYXRpb24gaXMgZXhwbGFpbmVkIGJ5IHRoZSBwcmVkaWN0aXZlIGF4ZXMgb2YgUExTIG9yIHRoZSBjb25zdHJhaW5lZCBwb3J0aW9uIG9mIFJEQSB0byBpbnRlcnByZXQgdGhlc2UgcmVzdWx0cy4KCiMgMy4gUmVkIEhlcnJpbmc/IChBS0EgTmVlZGxlIGluIGEgaGlnaGx5IG11bHRpY29sbGluZWFyIGhheXN0YWNrLikKCk9rLCB0aGlzIGlzIHJpZGljdWxvdXMsIGJ1dCBoZXJlIGlzIGEgZGF0YXNldCB3aXRoIDMwLCBzdHJvbmdseSBtdWx0aWNvbGxpbmVhciB2YXJpYWJsZXMgdGhhdCBkb24ndCBkaXNjcmltaW5hdGUgYmV0d2VlbiBncm91cHMgYW5kIG9ubHkgNSB2YXJpYWJsZXMgdGhhdCBkbyBkaXNjcmltaW5hdGUgYmV0d2VlbiBncm91cHMuCgpgYGB7cn0KZGYzIDwtIG1hcCgxOm5wZXJtLAogICAgICAgICAgICAgfmJhc2VfZGYgJT4lIAogICAgICAgICAgICAgICBzaW1fY292YXIoMzAsIDEsIDAuNywgbmFtZSA9ICJjb3JyIikgJT4lIAogICAgICAgICAgICAgICBzaW1fZGlzY3IoNSwgMSwgMCwgbmFtZSA9ICJkaXNjciIsIGdyb3VwX21lYW5zID0gYygtMC41LCAwLjUpKQopCgpkZjNbWzFdXSAlPiUgCiAgc2VsZWN0KC1ncm91cCkgJT4lIAogIGNvcigpICU+JSAKICBpaGVhdG1hcChyb3dfbGFiZWxzID0gVFJVRSwgY29sX2xhYmVscyA9IFRSVUUsIG5hbWUgPSAiY29yIikKYGBgCgoqKlBDQSAzKioKCmBgYHtyfQpwdmFscy5wY2EzIDwtIG1hcF9kYmwoZGYzLCBQQ0EudC50ZXN0KQptZWFuKHB2YWxzLnBjYTMpCmBgYAoKKipSREEgMioqCgpgYGB7cn0KcHZhbHMucmRhMyA8LSBtYXBfZGJsKGRmMywgUkRBLmFub3ZhKQptZWFuKHB2YWxzLnJkYTMpCmBgYAoKKipQTFMgMioqCgpgYGB7cn0Kb3V0cHV0LnBsczMgPC0gbWFwKGRmMywgUExTLnRlc3QpICU+JSBiaW5kX3Jvd3MoKQpgYGAKCiMjIFJlc3VsdHMgMwoKYGBge3J9CnJlc3VsdHMzIDwtIG91dHB1dC5wbHMzICU+JSBhZGRfY29sdW1uKHB2YWxzLnJkYTMsIHB2YWxzLnBjYTMpCnJlc3VsdHMzICU+JSBzdW1tYXJpemVfYWxsKGZ1bnMobWVhbiA9IG1lYW4pKQptZWFuUSA8LSBtZWFuKHJlc3VsdHMzJGBRMihjdW0pYCkKYGBgCgpgYGB7cn0KcHZhbHMzIDwtIHJlc3VsdHMzICU+JSAKICBzZWxlY3QoYFBMUyBwUTJgID0gcFEyLCBgUkRBIEFOT1ZBYCA9IHB2YWxzLnJkYTMsIGBQQ0EgdC10ZXN0YCA9IHB2YWxzLnBjYTMpICU+JSAKICBnYXRoZXIoa2V5ID0gIk1ldGhvZCIsIHZhbHVlID0gInAgdmFsdWUiKQoKZ2dwbG90KHB2YWxzMywgYWVzKHggPSBgcCB2YWx1ZWAsIGZpbGwgPSBNZXRob2QpKSArCiAgIyBnZW9tX2RlbnNpdHkoYWxwaGEgPSAwLjMpICsgCiAgZ2VvbV9oaXN0b2dyYW0ocG9zaXRpb24gPSAiZG9kZ2UiLCBiaW5zID0gMjApICsKICBnZW9tX3ZsaW5lKGFlcyh4aW50ZXJjZXB0ID0gMC4wNSksIGNvbG9yID0gInJlZCIpCmBgYAoKRXZlbiBpbiB0aGlzIHJpZGljdWxvdXMgc2NlbmFyaW8sIFBMUyBpcyBhYmxlIHRvICJmaW5kIiB0aGUgNSBkaXNjcmltaW5hdGluZyB2YXJpYWJsZXMgbW9zdCBvZiB0aGUgdGltZSBhbmQgZ2l2ZSBhIHNpZ25pZmljYW50IHAtdmFsdWUuICBUaGUgbWVhbiAkUV4yJCB2YWx1ZSBpcyBgciByb3VuZChtZWFuUSwgMilgLCBzbyBpdCdzIHBvc3NpYmxlIHRoYXQgc29tZSBvZiB0aG9zZSAic2lnbmlmaWNhbnQiIHAtdmFsdWVzIHdlcmUgb2J0YWluZWQgZnJvbSBtb2RlbHMgd2l0aCBwb29yIHByZWRpY3RpdmUgcG93ZXIuICBCZWNhdXNlIHdpdGhpbiBncm91cCBjby12YXJpYXRpb24gaXMgaGlnaCBhbmQgYmV0d2Vlbi1ncm91cCBjby12YXJpYXRpb24gaXMgbG93LCB0aGUgRi10ZXN0IHVzZWQgdG8gYW5hbHl6ZSBSREEgcmVzdWx0cyBnaXZlcyBhIGhpZ2ggcC12YWx1ZS4gIEkgdGhpbmsuIAoKCg==